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TECHNICAL  REPORT 
for  the  period  June  1,  1991— September  30,  1994 
#  AFOSR-91-0233 

Mathematical  Modeling  in  Plasma  Physics 
Natalia  Sternberg 

Department  of  Mathematics  and  Computer  Science 
Clark  University,  Worcester,  MA  01610 


The  proposed  project  consisted  of  three  parts: 

I.  To  study  the  transition  from  the  constant  ion  mobility  regime  at  the  center  of  the 
plasma  to  the  space-charge  regime  at  the  plasma  boundary  in  the  static  bounded  plasma 
problem. 

II.  To  understand  and  describe  both  analytically  and  numerically  the  characteristics 
for  the  dynamic  nonlinear  large  amplitude  perturbation  model  of  the  plasma  sheath  that 
has  been  developed  by  V.  Godyak  as  the  basis  of  the  theory  of  radio  frequency  discharges. 

III.  To  find  a  mathematical  model  for  determining  the  small  signal  capacitance  of  the 
plasma  sheath  at  intermediate  frequencies. 

A  number  of  problems  for  a  dynamic  model  of  the  sheath  have  been  successfully 
addressed  by  the  PI  in  collaboration  with  V.  Godyak  [1,2].  The  intent  of  Part  II  of  the 
proposal  was  to  modify  the  model  and  the  algorithm  in  [1,2]  so  it  becomes  appHcable  for 
Argon,  and  to  compare  the  numerical  and  experimental  data,  in  order  to  check  the  model 
and  the  theory  further.  During  the  period  covered  in  this  report,  the  goals  of  this  part  of 
the  project  have  been  a.chjeved: 

1.  The  investigator  was  able  to  give  a  complete  mathematical  analysis  for  a  self- 
consistent  dynamic  model  for  rf  sheaths  in  the  frequency  range  between  the  ion  and  electron 
plasma  frequencies  for  arbitrary  collision  parameters  and  arbitrary  rf  sheath  voltages. 

2.  Based  on  the  above  analysis,  the  investigator  developed  an  algorithm  which  permits 
a  complete  numerical  solution  of  the  problem  in  its  full  generality  without  any  restrictive 
assumptions. 

3.  Using  the  developed  algorithm,  the  investigator  has  computed  the  sheath  charac¬ 
teristics  for  Argon  and  Helium. 

4.  The  computed  results  were  compared  with  those  found  in  known  rf  sheath  theories 
[4-8],  and  with  those  obtained  experimentally  at  GTE  Laboratories  for  a  wide  range  of 
discharge  conditions.  It  was  found  that  our  mathematical  model  matches  the  experiment 
better  than  any  of  the  existing  models.  It  is  the  most  complete  and  generalized  model  for 
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the  rf  sheath  in  capacitive  rf  discharges. 


Based  on  the  obtained  results,  two  papers  have  been  written  by  the  PI: 

“A  comparison  of  rf  electrode  sheath  models”  (with  V.  Godyak  and  R.  Piejak),  IEEE 
Trans,  of  Plasma  Sci.,  21,  4,  1993 

“Solving  the  mathematical  model  of  the  electrode  sheath  in  symmetrically  driven  rf 
discharges”  (with  V.  Godyak),  J.  Comp.  Phys.,  Ill,  2,  347-353,  1994. 

In  the  first  paper  analytical  expressions  for  the  electrode  sheath  width  were  obtained 
for  different  time-averaged  and  dynamic  models  of  capacitive  rf  discharges.  It  has  been 
found  that  in  all  models  under  consideration,  the  sheath  widths  practically  coincide  if  they 
are  expressed  in  terms  of  the  dc  sheath  voltage  and  the  electrode  ion  current  to  the  rf 
electrode. 

In  the  second  paper,  a  complete  mathematical  analysis  was  given  for  a  self-consistent 
dynamic  model  of  rf  sheaths  in  the  frequency  range  between  the  plasma  ion  and  elec¬ 
tron  frequency.  Arbitrary  collision  parameters  and  arbitrary  rf  sheath  voltages  have  been 
considered.  Based  on  this  analysis  an  algorithm  has  been  developed  to  solve  numerically 
the  mathematical  problem  in  its  full  generality.  The  developed  method  was  then  used  to 
obtain  the  electrode  sheath  characteristics  and  to  study  their  dependence  on  the  rf  sheath 
conductivity  parameter  for  Argon. 

In  somewhat  more  detail,  the  results  obtained  to  date  can  be  summarized  as  follov.-s. 

A  self-consistent  hydrodynamic  model  for  the  sheaths  of  a  symmetrically  driven  rf 
discharge  was  formulated  in  Ref.  [3].  In  Ref.  [1],  this  model  was  generalized.  In  the 
model,  under  the  action  of  the  rf  field,  the  plasma-sheath  interface  exhibits  an  oscillatory 
behavior.  Therefore,  the  position  of  the  plasma-sheath  interface  at  any  given  time  6  can  be 
described  by  a  periodic  function  X{6)  that  achieves  its  maximum  Aj  at  some  time  d  —  di 
and  its  minimum  A2  at  6  =  ^2-  To  the  right  of  the  moving  interface,  in  the  plasma, 
the  electron  density  is  equal  to  the  ion  density  (plasma  neutrality).  To  the  left  of  the 
moving  interface,  in  the  sheath,  the  electron  density  is  equal  to  zero  (space  charge).  From 
the  physical  point  of  view  the  generalized  model  [1]  is  the  most  accurate  and  justified, 
mathematically,  however,  it  leads  to  a  complicated  moving  boundary  problem  that  in  the 
past  has  had  to  be  simplified  in  order  to  be  solved.  A  number  of  such  simplified  models 
appeared  in  the  literature  (see  [4-8]),  but  the  mathematical  simplification  have  lead  to 
physical  inconsistencies  and  inaccuracies.  In  Ref.  [1],  we  worked  out  the  generalized 
model  under  the  single  restrictive  assumption  that  the  function  A  is  symmetric  with  respect 
to  its  minimum  position  A2.  A  comparison  of  the  analytical  and  numerical  solutions  of 
the  mathematical  model  with  the  available  experimental  results  for  mercury  vapor  have 
demonstrated  a  good  agreement  in  corresponding  values.  The  mathematical  analysis, 
which  underlies  the  numerical  results,  can  be  found  in  [2].  In  general,  however,  A  is 
not  symmetric  with  respect  to  its  minimum  position  A2,  and  always  imposing  such  an 
assumption  might  lead  to  inaccuracies  for  other  gases.  Indeed,  such  an  assumption  leads 
to  corresponding  symmetries  in  the  solutions  to  the  problem.  Such  symmetry,  however,  is 
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not  to  be  expected.  In  fact,  asymmetry  becomes  more  pronounced  for  large  oscillations  of 
the  plasma-sheath  interface.  How  to  solve  the  moving  boundary  problem,  which  describes 
the  generalized  model  [1],  without  any  restrictive  assumptions  or  simplification  has  up 
until  now  remained  an  open  question.  We  were  able  to  answer  this  question  by  giving 
a  detailed  mathematical  analysis  of  the  generalized  model.  This  enabled  us  to  solve  the 
problem  numerically  in  its  full  generality,  and  to  obtain  the  characteristics  of  the  electrode 
sheath.  Using  the  numerical  results,  we  were  also  able  to  obtain  an  approximation  formula 
for  the  sheath  characteristics  for  arbitrary  sheath  voltages  and  collision  parameters.  This 
formula  is  similar  to  the  Child-Langmuir  Law  and  can  easily  be  used  for  the  experimental 
modeling  of  rf  discharges. 

In  most  applications,  one  has  to  deal  with  bounded  plasmas.  Our  approach  can  be 
used  for  any  bounded  pla.sma  created  by  rf  antennas  or  rf  electrodes. 

The  above  results  led  to  the  following  new  problem  that  originally  has  not  been  pro¬ 
posed.  In  order  to  solve  the  model  for  the  sheath,  one  has  to  specify  the  boundary  condition 
at  the  plasma-sheath  interface  which  would  provide  a  smooth  transition  from  the  plasma  to 
the  sheath.  Since  there  is  no  strict  boundary  between  the  plasma  and  the  sheath,  it  is  not 
obvious  how  to  determine  such  a  boundary  condition.  Most  frequently  the  electric  field  at 
the  plasma-sheath  boundary  Ei  is  chosen  to  be  zero.  This  choice  results  in  discontinuities 
of  the  gradients  of  the  potential,  plasma-density,  and  velocity,  and  in  some  cases  it  leads 
to  an  infinite  sheath  width.  Several  attempts  have  been  made  to  remove  these  disconti¬ 
nuities  by  introducing  transition  layers.  Another  way  to  chose  the  boundary  condition  at 
the  plasma-sheath  interface  is  to  specify  the  value  of  the  electric  field  where  the  plasma 
neutrality  is  violated  by  setting  Ei  =  T^/Xoii  where  Tg  is  the  electron  temperature  and 

is  the  electron  Debye  radius  at  the  plasma-sheath  boundary.  This  nonzero  boundary 
condition  for  the  electric  field  has  been  introduced  by  Godyak  [9],  and  has  previously  been 
used  by  the  PI  to  solve  the  static  sheath  model  and  the  dynamic  sheath  model.  However, 
the  precise  value  for  the  electric  field  at  the  plasma-sheath  interface  is  still  being  disputed. 
This  dispute  leads  to  the  question  about  the  effects  of  choice  of  boundary  conditions  on 
the  sheath  characteristics.  This  question  has  been  addressed  by  the  PI  in  the  paper 

“Effects  of  boundary  conditions  in  the  dynamic  model  of  the  rf  sheath”  (with  V. 
Godyak),  IEEE  Trans.  Magnetics,  30,  5,  1994. 

Some  of  the  results  in  this  paper  were  presented  at  the  Conference  on  the  Compu¬ 
tations  of  Electromagnetic  Fields  (COMPUMAG)  in  Miami  on  November  5-6,  1993,  and 
appeared  in  the  COMPUMAG  Conference  Records,  1993. 

In  that  paper,  the  calculations  were  performed  for  the  benchmark  modeling  case  of 
the  symmetrically  driven  rf  discharge  in  Argon  at  different  values  A  of  the  normalized 
electric  field  at  the  plasma-sheath  interface;  A  =  0;  —0.1;  —0.3;  —0.5;  —0.7;  —0.9;  — 1. 
The  results  that  have  been  obtained  show  that  even  for  relatively  high  rf  sheath  voltages, 
both  the  dc  sheath  voltage  and  the  rf  sheath  capacitance  depend  in  a  significant  way  on 
the  choice  of  the  value  for  the  electric  field  at  the  plasma-sheath  interface.  These  results 
are  surprising  since  for  high  rf  sheath  voltages,  it  has  always  been  believed  that  the  sheath 
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the  sheath  characteristics  behave  asymptotically  and  are  not  sensitive  to  the  choice  of  the 
boundary  conditions.  Thus,  the  question  about  the  choice  of  boundary  condition  at  the 
plasma-sheath  interface  for  an  adequate  rd  sheath  model  is  indeed  important  and  relevant. 
An  analysis  given  by  Riemann  in  [10]  suggests  that  \A\  ^  (Aj^/T),  where  L  is  the  plasma 
width.  For  rf  discharges  encountered  in  applications,  this  value  is  much  smaller  than  1.  On 
the  other  hands,  Godyak’s  boundary  condition  suggests  that  A  =  —  1.  It  remained  to  be 
investigated  which  values  of  A  are  appropriate.  Such  an  investigation  has  been  conducted 
by  the  PI  for  the  static  boimded  plasma  problem  (see  Part  I).  The  obtained  results  will 
be  described  below. 

During  the  period  covered  by  this  report,  the  PI  has  also  been  working  on  the  static 
bounded  plasma  problem  with  finite  nonneutrality.  A  number  of  problems  concerning 
the  static  bounded  plasma  problem  have  been  addressed  by  the  PI  in  collaboration  with 
Godyak  in  [11].  In  that  paper,  the  hydrodynamic  models  for  the  plasma  and  the  sheath 
have  been  studied  separately  using  Godyak’s  boundary  condition.  In  the  plasma  model, 
one  assumes  that  the  ion  density  equals  the  electron  density  (qucisineutrality).  In  actual 
applications,  this  assumption  can  be  made  with  sufficiently  high  accuracy  in  the  center 
of  the  plasma.  Closer  to  the  boundary  this  is  no  longer  true.  In  the  sheath  model,  one 
assumes  that  the  electron  density  equals  zero.  This  assumption  can  be  justified  only  very 
close  to  the  wall.  These  facts  immediately  raise  an  important  questions.  In  what  region  can 
the  solution  of  the  total  bounded  plasma  problem  be  approximated  by  the  solutions  of  the 
plasma  and  the  sheath  models?  In  order  to  answer  this  question,  a  mathematical  analysis 
of  the  static  model  has  been  given,  and  all  the  needed  algorithms  have  been  developed. 
The  following  results  for  the  collisionless  case  with  the  nonneutrality  parameter  q  =  0.01 
have  been  obtained; 

1.  Using  an  iteration  method  we  have  found  the  position  of  the  wall:  =  0.6558. 

2.  The  plasma  model  approximates  the  bounded  plasma  problem  in  the  region 
a:  <  0.56. 

3.  The  sheath  model  approximates  the  bounded  plasma  problem  in  the  region 
0.63  <  X  <  Xyj, 

4.  In  the  region  0.56  <  a:  <  0.63,  the  bounded  plasma  problem  cannot  be  approxi¬ 
mated  by  either  the  plasma  or  the  sheath  model.  This  region  corresponds  to  the  presheath 
region. 

These  results  lead  to  the  following  important  question.  What  position  of  the  plasma- 
sheath  interface  is  specified  by  Godyak’s  boundary  condition,  and  what  position  is  specified 
by  Riemann’s  boundary  condition?  In  order  to  study  this  problem,  we  have  solved  the 
plasma  and  the  sheath  models  varying  the  boundary  condition  at  the  plasma-sheath  in¬ 
terface.  Comparing  the  so  obtained  solutions  with  the  solution  of  the  bounded  plasma 
problem  the  following  has  been  found: 

5.  Godyak’s  boundary  condition  specifies  the  position  at  the  plasma-sheath  interface 
at  X  =  0.63.  Therefore,  this  boundary  condition  can  indeed  be  used  as  an  initial  con- 
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dition  for  the  sheath  model  to  approximate  the  bounded  plasma  problem  in  the  region 
0.63  <  X  <  Xtjj. 

6.  Riemann’s  boundary  condition  specifies  the  position  at  the  plasma-sheath  interface 
at  a;  =  0.56.  This  boundary  condition  can  be  used  as  a  boundary  condition  for  the  plasma 
model  to  approximate  the  bounded  plasma  problem  in  the  region  0  <  a;  <  0.56. 

It  becomes  obvious  from  the  above  results  that  one  caimot  just  combine  the  solutions 
of  the  plasma  and  the  sheath  models  to  obtain  an  accurate  approximation  of  the  solution  of 
the  bounded  plasma  problem.  In  order  to  find  an  accurate  approximation  of  the  bounded 
plasma  problem  one  would  have  to  model  not  only  the  plasma  and  the  sheath,  but  also 
the  presheath.  This  conclusion  lead  to  another  question.  How  relevant  is  the  presheath 
in  applications?  In  order  to  answer  this  question,  we  have  computed  the  low  frequency 
capacitance  using  the  total  plasma  model,  and  using  the  sheath  model.  A  comparison  of 
the  obtained  results  has  shown  the  following; 

7.  In  some  applications  the  presheath  can  be  ignored  and  the  bounded  plasma  problem 
can  indeed  be  modeled  by  the  plasma  and  the  sheath  models  separately.  The  so  obtained 
integral  characteristics  of  the  problem  are  within  the  accuracy  of  the  experiment. 

The  preliminary  findings  concerning  the  static  bounded  plasma  problem  were  pre¬ 
sented  at  the  IEEE  International  Conference  on  Plasma  Science,  Santa  Fe,  New  Mexico, 
June,  1994.  A  detailed  description  of  the  above  results  is  currently  being  prepared  for 
publication. 


During  the  period  covered  by  this  report,  the  following  papers  have  been  published 
by  the  PI: 

1.  Rf  sheath  parameters  for  Argon  and  Helium  discharges.  Conference  Records, 
IEEE  International  Conference  on  Plasma  Science,  Williamsburg,  VA,  June,  1991  (with 
V.  Godyak  and  R.  Piejak). 

2.  Sheath  capacitance  at  low  and  high  frequencies  (with  V.  Godyak),  Proceed.  XX 
Int.  Conf.  on  Phenomena  in  Ionized  Gases,  Pisa,  Italy,  July,  1991. 

3.  A  Hartman-Grobman  Theorem  for  maps.  Proceed.  Int.  Conf.  on  Theory  and  Appl. 
of  Diff.  Eq.,  Edinburg,  Texas,  May,  1991,  in  “Ordinary  and  delay  differential  equations”, 
Eds.  J.  Hale  and  J.  Wiener,  Pitman  Research  Notes  in  Mathematics  Series,  1992. 

4.  A  moving  boundary  value  problem  in  Plasma  Physics  (with  V.  Godyak),  Proceed. 
Ill  Potsdam-V  Kiev  Int.  Workshop  on  Nonlinear  Processes  in  Physics,  Potsdam,  NY,  Au¬ 
gust,  1991,  in  “Nonlinear  Processes  in  Physics”,  Eds.  A.S.  Fokas,  D.J.Kaup,  A. C.  Newell, 
V.E.  Zakharov,  Springer- Verlag,  1993. 

5.  A  Hartman-Grobman  Theorem  for  a  class  of  retarded  functional  differential  equa¬ 
tions,  J.  Math.  Analysis  and  Appl.,  175,  1993. 


5 


6.  A  comparison  of  rf  electrode  sheath  models  (with  V.  Godyak  and  R.  Piejak),  IEEE 
Trans,  of  Plasma  Sci.,  21,  4,  1993. 

7.  Solving  the  mathematical  model  of  the  electrode  sheath  in  symmetrically  driven  rf 
discharges  (with  V. A.  Godyak),  J.  Comp.  Phys.,  Ill,  2,  347-353,  1994. 

8.  Effects  of  boundary  conditions  in  the  dynamic  model  of  the  rf  sheath,  (with 
V.A.  Godyak),  IEEE  Trans.  Magnetics,  30,  5,  3100-3103,  1994. 

During  the  period  covered  by  this  report,  the  following  presentations  at  international 
conferences  have  been  given  by  the  PI: 

1.  Rf  sheath  parameters  for  Argon  and  HeUum  discharges,  IEEE  International  Con¬ 
ference  on  Plasma  Science,  Williamsbmrg,  VA,  June,  1991. 

2.  Sheath  capacitance  at  low  and  high  frequencies  (with  V.  Godyak),  Proceed.  XX 
Int.  Conf.  on  Phenomena  in  Ionized  Gases,  Pisa,  Italy,  Jiily,  1991. 

3.  Mathematical  modeling  in  Plasma  Physics,  invited  participant.  III  Potsdam-V  Kiev 
International  Workshop  on  Nonlinear  Processes  in  Physics,  Clarkson  University,  Potsdam, 

New  York,  August,  1991. 

4.  A  Hartman-Grobman  theorem  for  maps,  SIAM  Conference  on  Applications  of 
Dynamical  Systems,  Snowbird,  Utah,  October,  1992. 

5.  Effects  of  boundary  conditions  in  the  dynamic  model  of  the  rf  sheath,  COMPUMAG 
9th  Conference  on  the  Computation  of  Electromagnetic  Fields,  Miami,  Florida,  November, 
1993. 

6.  Approximation  of  the  bounded  plasma  model  by  the  plasma  and  the  sheath  models, 
IEEE  Int.  Conf.  on  Plasma  Science,  Santa  Fe,  New  Mexico,  June,  1994. 
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